#include "sweSolver_valuesTestFunctions.hpp"
#include "settingsFunctions.hpp"
#include "meshTypeFunctions.hpp"
#include "bc_base_all.hpp"
#include "bc_swe_all.hpp"

#include "streamAccess.hpp"

#include "catch2_testDir.hpp"
#include <catch2/catch_test_macros.hpp>

using namespace hms;
namespace fs = std::filesystem;

template<typename MeshType>
MeshType vShapedDomain_mesh(){
	// std::shared_ptr<spdlog::logger> logger { hms::log::getLogger() };
	// logger->set_level(hms::log::Level::trace);

	constexpr Index
		nx {8},
		ny {8},
		ng {1};
	/* a v-shape is easier to achieve with the center of the mesh being the
	 * coordinate origin */
	constexpr scalar
		dx {1},
		dy {1},
		offsetX {-dx * 0.5 * nx},
		offsetY {-dy * 0.5 * ny};

	return uniMeshAs<MeshType>(nx, ny, dx, dy, true, offsetX, offsetY);
}

template<typename MeshType>
void vShapedDomain_settings( SWESolver<MeshType>& solver, Time& time ){

	// time.writeTimeInterval(1);
	time.endAt( std::numeric_limits<scalar>::max() );
	time.writeTimeInterval( std::numeric_limits<scalar>::max() );
	time.maxSteps(2);

	Field<MeshType>
		& d { solver.d() },
		& v { solver.v() },
		& z { solver.z() };
	
	v.values().colwise() = Array2s{0.5, 0};
	d.values() = 1;
	z.values() = 0;

	/* create a v-shape, with the crease being a horizontal line through the
	 * center of the mesh */
	constexpr scalar
		slopeX {-0.1},
		slopeY {0.1};

	const MeshType& mesh { solver.mesh() };

	mesh.forEachCell([&](const auto& cell){
		Index i { cell.index() };
		Vector2s centroid { cell.centroid() };
		/* with the next line, we get an inclined plane */
		scalar zVal { centroid.x() * slopeX };
		/* now add the v-shape */
		zVal += std::abs( centroid.y() ) * slopeY;
		solver.z().values()(0,i) = zVal;
	});

	const BoundaryPatch
		& bottom { mesh.boundaryPatch("bottom") },
		& top    { mesh.boundaryPatch("top") },
		& right  { mesh.boundaryPatch("right") },
		& left   { mesh.boundaryPatch("left") };

	scalar vGhost {1};
	d.setBoundaryCondition( OrthogonalFlowVelocity_waterdepth{d, left, v, vGhost} );
	v.setBoundaryCondition( OrthogonalFlowVelocity_velocity  {v, left,    vGhost} );

	d.setBoundaryCondition( FreeOutflow_waterdepth{d, right, false} );
	v.setBoundaryCondition( FreeOutflow_velocity  {v, right} );

	d.setBoundaryCondition( ZeroGradient {d, top} );
	v.setBoundaryCondition( Slip_velocity{v, top} );

	d.setBoundaryCondition( ZeroGradient {d, bottom} );
	v.setBoundaryCondition( Slip_velocity{v, bottom} );

	z.setBoundaryCondition( StaticValue{z, left  , false} );
	z.setBoundaryCondition( StaticValue{z, right , false} );
	z.setBoundaryCondition( StaticValue{z, bottom, false} );
	z.setBoundaryCondition( StaticValue{z, top   , false} );

	solver.settings().useFriction = true;
	solver.settings().edgeTraversal = EdgeTraversal::ordered_secondOrder;
	/* mandate second order reconstruction */
	solver.settings().useMuscl = true;
	solver.settings().limiterFunction = LimiterFunction::minMod;

	solver.initialise();
}

template<typename MeshType>
void vShapedDomain_checkResults( const SWESolver<MeshType>& solver ){

	const MeshType& mesh { solver.mesh() };

	Field<MeshType>
		dTargetField {mesh, 1},
		vTargetField {mesh, 2};

	ArrayXXs
		& dTarget { dTargetField.values() },
		& vTarget { vTargetField.values() };

	dTarget <<
0      , 1.04584, 1       , 1       , 1       , 1       , 1       , 1       , 1       , 0,
1.12144, 1.05942, 0.993935, 0.983084, 0.983084, 0.983084, 0.983084, 0.983084, 0.983084, 1,
1.12144, 1.07012, 1.00474 , 0.993892, 0.993892, 0.993892, 0.993892, 0.993892, 0.993892, 1,
1.10925, 1.06564, 1.00063 , 0.989781, 0.989781, 0.989781, 0.989781, 0.989781, 0.989781, 0.989153,
1.13363, 1.11017, 1.04409 , 1.03324 , 1.03324 , 1.03324 , 1.03324 , 1.03324 , 1.03324 , 1.01085 ,
1.13363, 1.11017, 1.04409 , 1.03324 , 1.03324 , 1.03324 , 1.03324 , 1.03324 , 1.03324 , 1.01085 ,
1.10925, 1.06564, 1.00063 , 0.989781, 0.989781, 0.989781, 0.989781, 0.989781, 0.989781, 0.989153,
1.12144, 1.07012, 1.00474 , 0.993892, 0.993892, 0.993892, 0.993892, 0.993892, 0.993892, 1,
1.12144, 1.05942, 0.993935, 0.983084, 0.983084, 0.983084, 0.983084, 0.983084, 0.983084, 1,
0      , 1.04584, 1       , 1       , 1       , 1       , 1       , 1       , 1       , 0;
	vTarget.reshaped() <<
0.5, 0, 0.772498, -0.129072,  0.634721, -0.134989,  0.634721, -0.134989,  0.634721, -0.134989,  0.634721, -0.134989,  0.634721, -0.134989,  0.634721, -0.134989,  0.634721, -0.134989,  0.5,        0,
1  , 0, 0.986996,  0.242423,  0.793723,  0.258266,  0.759337,  0.259905,  0.759337,  0.259905,  0.759337,  0.259905,  0.759337,  0.259905,  0.759337,  0.259905,  0.759337,  0.259905,  0.634721,   0.134989,
1  , 0, 0.984842,  0.259641,  0.792006,  0.275711,  0.757975,  0.277523,  0.757975,  0.277523,  0.757975,  0.277523,  0.757975,  0.277523,  0.757975,  0.277523,  0.757975,  0.277523,  0.634721,   0.134989,
1  , 0, 0.985376,  0.255723,  0.792467,  0.273828,  0.758576,  0.27532,   0.758576,  0.27532,   0.758576,  0.27532,   0.758576,  0.27532,   0.758576,  0.27532,   0.758576,  0.27532,   0.636203,   0.169525,
1  , 0, 0.977511,  0.0726889, 0.786238,  0.0772342, 0.753177,  0.0777655, 0.753177,  0.0777655, 0.753177,  0.0777655, 0.753177,  0.0777655, 0.753177,  0.0777655, 0.753177,  0.0777655, 0.633271,   0.0327553,
1  , 0, 0.977511, -0.0726889, 0.786238, -0.0772342, 0.753177, -0.0777655, 0.753177, -0.0777655, 0.753177, -0.0777655, 0.753177, -0.0777655, 0.753177, -0.0777655, 0.753177, -0.0777655, 0.633271,  -0.0327553,
1  , 0, 0.985376, -0.255723,  0.792467, -0.273828,  0.758576, -0.27532,   0.758576, -0.27532,   0.758576, -0.27532,   0.758576, -0.27532,   0.758576, -0.27532,   0.758576, -0.27532,   0.636203,  -0.169525,
1  , 0, 0.984842, -0.259641,  0.792006, -0.275711,  0.757975, -0.277523,  0.757975, -0.277523,  0.757975, -0.277523,  0.757975, -0.277523,  0.757975, -0.277523,  0.757975, -0.277523,  0.634721,  -0.134989,
1  , 0, 0.986996, -0.242423,  0.793723, -0.258266,  0.759337, -0.259905,  0.759337, -0.259905,  0.759337, -0.259905,  0.759337, -0.259905,  0.759337, -0.259905,  0.759337, -0.259905,  0.634721,  -0.134989,
0.5, 0, 0.772498,  0.129072,  0.634721,  0.134989,  0.634721,  0.134989,  0.634721,  0.134989,  0.634721,  0.134989,  0.634721,  0.134989,  0.634721,  0.134989,  0.634721,  0.134989,  0.5,        0;


	scalar prec {1e-5};
	const ArrayXXs
		& d { solver.d().values() },
		& v { solver.v().values() };
	bool
		dDomainInRange { solver.d().domainValues()
			.isApprox(dTargetField.domainValues(), prec)
		},
		vDomainInRange { solver.v().domainValues()
			.isApprox(vTargetField.domainValues(), prec)
		},
		dInRange { d.isApprox(dTarget, prec) },
		vInRange { v.isApprox(vTarget, prec) };

	Index
		ng  { mesh.nGhostLayers() },
		nxg { 2*ng + mesh.nCellsX() },
		nyg { 2*ng + mesh.nCellsY() };

	auto reshape = [&](const auto& arr){
		return arr.reshaped( nxg * arr.rows(), nyg );
	};

	if ( !dDomainInRange || !vDomainInRange ){
		std::cout << "Domain values not in range. Printing...\n";

		auto printIfNotEqual = [&](
			const auto& target, const auto& actual, Index cell,  char field
		){
			if ( !target.col(cell).isApprox( actual.col(cell) ) )
				std::cout
					<< "Cell #" << cell
					<< ", field " << field
					<< "\ntarget: " << target.col(cell).transpose()
					<< "\nactual: " << actual.col(cell).transpose()
					<< '\n';
		};

		forEachCell<WhichCells::domain, true, true>(
			mesh, [&](Index cell){
				printIfNotEqual( dTarget, d, cell, 'd' );
				printIfNotEqual( vTarget, v, cell, 'v' );
			}
		);
	}
	REQUIRE( dDomainInRange );
	REQUIRE( vDomainInRange );

	if (!dInRange || !vInRange){
		std::cout << "Boundary values not in range. Printing...\n";

		auto printBoundaryIfNotEqual = [&](
			const Field<MeshType>& target,
			const Field<MeshType>& actual,
			const BoundaryPatch& p,
			char field
		){
			const auto
				& tVals { target.ghostValues(p) },
				& aVals { actual.ghostValues(p) };
			if ( !tVals.isApprox(aVals, prec) ){
				std::cout
					<< "Field \"" << field << "\", patch \"" << p.name()
					<< "\", error " << (tVals - aVals).matrix().norm()
					<< " > " << prec
					<< "\ntarget\tactual\n"
					<< concat_horizontal( tVals.reshaped(), aVals.reshaped() )
					<< '\n';
			}
		};

		for ( const BoundaryPatch& p : mesh.boundaryPatches() ){
			printBoundaryIfNotEqual( dTargetField, solver.d(), p, 'd' );
			printBoundaryIfNotEqual( vTargetField, solver.v(), p, 'v' );
		}
	}

	REQUIRE( dInRange );
	REQUIRE( vInRange );
}


template<typename MeshType>
MeshType dambreakTest_mesh(){

	constexpr Index
		nx {30},
		ny {4},
		ng {1};
	/* A discontinuity in the water depth is later created at center of the mesh.
	 * Therefore, it is helpful to have the center be the coordinate origin. */
	constexpr scalar
		dx {0.2},
		dy {0.2},
		offsetX {-dx * 0.5 * nx},
		offsetY {-dy * 0.5 * ny};

	return uniMeshAs<MeshType>(nx, ny, dx, dy, true, offsetX, offsetY);
}

template<typename MeshType>
void dambreakTest_settings(
	SWESolver<MeshType>& solver,
	Time& time,
	bool useFriction
){

	// time.writeTimeInterval(1);
	// time.endAt( std::numeric_limits<scalar>::max() );
	time.endAt(0.1);
	time.writeTimeInterval( std::numeric_limits<scalar>::max() );
	// time.maxSteps(4);

	Field<MeshType>
		& d { solver.d() },
		& v { solver.v() },
		& z { solver.z() };
	
	v.values() = 0;
	d.values() = 4;
	z.values() = 0;

	const MeshType& mesh { solver.mesh() };

	/* all cells in the right half of the mesh are set dry */
	solver.mesh().forEachCell( [&](const auto& cell){
		if ( cell.centroid().x() > 0 )
			d.values()( 0, cell.index() ) = 0;
	} );

	const BoundaryPatch
		& bottom { mesh.boundaryPatch("bottom") },
		& top    { mesh.boundaryPatch("top") },
		& right  { mesh.boundaryPatch("right") },
		& left   { mesh.boundaryPatch("left") };

	d.setBoundaryCondition( FreeOutflow_waterdepth{d, left, false} );
	v.setBoundaryCondition( FreeOutflow_velocity  {v, left} );

	d.setBoundaryCondition( FreeOutflow_waterdepth{d, right, false} );
	v.setBoundaryCondition( FreeOutflow_velocity  {v, right} );

	d.setBoundaryCondition( ZeroGradient {d, top} );
	v.setBoundaryCondition( Slip_velocity{v, top} );

	d.setBoundaryCondition( ZeroGradient {d, bottom} );
	v.setBoundaryCondition( Slip_velocity{v, bottom} );

	z.setBoundaryCondition( StaticValue{z, left  , false} );
	z.setBoundaryCondition( StaticValue{z, right , false} );
	z.setBoundaryCondition( StaticValue{z, bottom, false} );
	z.setBoundaryCondition( StaticValue{z, top   , false} );

	auto& settings { solver.settings() };

	settings.useFriction = useFriction;
	settings.edgeTraversal = EdgeTraversal::ordered_secondOrder;
	/* mandate second order reconstruction */
	settings.useMuscl = true;
	settings.limiterFunction = LimiterFunction::minMod;
	settings.maxCourant = 0.3;

	solver.initialise();
}

template<typename MeshType>
void dambreakTest_checkResults(
	const SWESolver<MeshType>& solver,
	bool useFriction
){

	const MeshType& mesh { solver.mesh() };
	Index
		nx { mesh.nCellsX() },
		ny { mesh.nCellsY() },
		ng { mesh.nGhostLayers() };

	fs::path resultsPath { globalVar::testDir };
	if (useFriction)
		resultsPath /= "sweSolver_valuesTest_dambreak_friction.txt";
	else
		resultsPath /= "sweSolver_valuesTest_dambreak_noFriction.txt";

	REQUIRE( fs::is_regular_file(resultsPath) );

	std::ifstream resultsFile { getReadFile(resultsPath) };


	REQUIRE( 'v' == read<char>(resultsFile) );
	Array2Xs vRow (2, nx);
	read(resultsFile, vRow);

	REQUIRE( 'd' == read<char>(resultsFile) );
	Array1Xs dRow (1, nx);
	read(resultsFile, dRow);

	auto getDomainBlock = [&](const Field<MeshType>& field){
		const Index& nc { field.nComponents() };
		return field.values()
			.reshaped(nc*(nx+2*ng), ny+2*ng)
			.block(nc*ng, ng, nc*nx, ny);
	};

	auto vDomain { getDomainBlock( solver.v() ) };
	auto dDomain { getDomainBlock( solver.d() ) };

	// std::cout
	// 	<< "\nvRow:\n" << vRow.row(0).head(8)
	// 	<< "\nvDomain, col 0:\n" << vDomain.col(0).reshaped(2,nx).row(0).head(8)
	// 	<< '\n';

	for ( Index i{0}; i<dDomain.cols(); ++i ){
		REQUIRE( vDomain.col(i).isApprox( vRow.reshaped() ) );
		REQUIRE( dDomain.col(i).isApprox( dRow.reshaped() ) );
	}
	// if constexpr ( std::is_same_v<MeshType, UniMesh> )
	// 	printAsStructMesh( solver.v().values(), solver.mesh() );
}


#define INSTANTIATE_FOR_MESH_TYPE(MeshType) \
template MeshType vShapedDomain_mesh(); \
template void vShapedDomain_settings( SWESolver<MeshType>&, Time& ); \
template void vShapedDomain_checkResults( const SWESolver<MeshType>& solver ); \
template MeshType dambreakTest_mesh(); \
template void dambreakTest_settings( SWESolver<MeshType>&, Time&, bool ); \
template void dambreakTest_checkResults( const SWESolver<MeshType>& solver, bool );

INSTANTIATE_FOR_MESH_TYPE( UniMesh)
INSTANTIATE_FOR_MESH_TYPE(RectMesh)

#undef INSTANTIATE_FOR_MESH_TYPE


// template  UniMesh vShapedDomain_mesh();
// template RectMesh vShapedDomain_mesh();

// template void vShapedDomain_settings( SWESolver< UniMesh>&, Time& );
// template void vShapedDomain_settings( SWESolver<RectMesh>&, Time& );

// template void vShapedDomain_checkResults( const SWESolver< UniMesh>& solver );
// template void vShapedDomain_checkResults( const SWESolver<RectMesh>& solver );